/***************************************************************************
 *   Copyright (C) 2008 by Eduardo Gurgel Pinho                            *
 *   edgurgel@gmail.com                                                    *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 3 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
***************************************************************************/

#include "quadrilateral.h"

#include "math/ray.h"
#include "material/object_material.h"

#include <cmath>
#include <limits>

namespace Sirat {

namespace Object {

Quadrilateral::Quadrilateral(Material::Object_Material* _object_material,const Math::Vector& _vertex0,const Math::Vector& _vertex1,const Math::Vector& _vertex2,const Math::Vector& _vertex3) : Object(_object_material) {
	vertices[0] = _vertex0;
	vertices[1] = _vertex1;
	vertices[2] = _vertex2;
	vertices[3] = _vertex3;
}

Quadrilateral::Quadrilateral() : Object() { }

Quadrilateral::Quadrilateral(const Quadrilateral& _quadrilateral) : Object(_quadrilateral) {
	vertices[0] = _quadrilateral.vertices[0];
	vertices[1] = _quadrilateral.vertices[1];
	vertices[2] = _quadrilateral.vertices[2];
	vertices[3] = _quadrilateral.vertices[3];
}

Quadrilateral::~Quadrilateral() { }

Math::Vector Quadrilateral::get_vertex(const unsigned short _index) const { return vertices[_index]; }

Math::Vector Quadrilateral::get_normal(const Math::Vector& _point) const { 
	const Math::Vector v1 = vertices[1] - vertices[0];
	const Math::Vector v2 = vertices[3] - vertices[0];
	return (v1.cross_product(v2)).normalize();
}

void Quadrilateral::set_vertex(const Math::Vector& _vertex, const unsigned short _index) { vertices[_index] = _vertex; }

tuple<double,Math::Vector,Object*> Quadrilateral::intersection(const Math::Ray& _ray) {
	//Ares Lagae and Philip Dutré ray-quadrilateral intersection

	
	std::numeric_limits<double> double_limit;
	const double infinity = double_limit.infinity();
	const double epsilon = double_limit.epsilon();
	//const double epsilon = 0.00001;

	Math::Vector intersection_point;
	double u,v;
	
	Math::Vector E_01 = vertices[1] - vertices[0];
	Math::Vector E_03 = vertices[3] - vertices[0];
	Math::Vector P = _ray.get_direction().cross_product(E_03);
	double det = E_01.dot_product(P);
	if (std::fabs(det) < epsilon) return make_tuple(infinity,intersection_point,this);
	double inv_det = 1.0/det;
	Math::Vector T = _ray.get_origin() - vertices[0];
	double alpha = T.dot_product(P) * inv_det;
	if (alpha < 0.0) return make_tuple(infinity,intersection_point,this);
	if (alpha > 1.0) return make_tuple(infinity,intersection_point,this);
	Math::Vector Q = T.cross_product(E_01);
	double beta = _ray.get_direction().dot_product(Q) * inv_det;
	if (beta < 0.0) return make_tuple(infinity,intersection_point,this);
	if (beta > 1.0) return make_tuple(infinity,intersection_point,this);

	if ((alpha + beta) > 1.0) {

    // Rejects rays that intersect the plane of Q either on the
    // left of the line V11V10 or on the right of the line V11V01.

		Math::Vector E_23 = vertices[3] - vertices[2];
		Math::Vector E_21 = vertices[1] - vertices[2];
		Math::Vector P_prime = _ray.get_direction().cross_product(E_21);
		double det_prime = E_23.dot_product(P_prime);
		if (std::abs(det_prime) < epsilon) return make_tuple(infinity,intersection_point,this);
		double inv_det_prime = 1.0 / det_prime;
		Math::Vector T_prime = _ray.get_origin() - vertices[2];
		double alpha_prime = T_prime.dot_product(P_prime) * inv_det_prime;
		if (alpha_prime < 0.0) return make_tuple(infinity,intersection_point,this);
		Math::Vector Q_prime = T_prime.cross_product(E_23);
		double beta_prime = _ray.get_direction().dot_product(Q_prime) * inv_det_prime;
		if (beta_prime < 0.0) return make_tuple(infinity,intersection_point,this);
	}
	// Compute the ray parameter of the intersection point, and
	// reject the ray if it does not hit Q.

	double t = E_03.dot_product(Q) * inv_det;
	if (t < 0.0) return make_tuple(infinity,intersection_point,this);
	else { 
		intersection_point = _ray.get_origin() + _ray.get_direction()*t;
		return make_tuple(t,intersection_point,this); 
	}
	// Compute the barycentric coordinates of the fourth vertex.
	// These do not depend on the ray, and can be precomputed
  	// and stored with the quadrilateral.  

	double alpha_11, beta_11;
	Math::Vector E_02 = vertices[2] - vertices[0];
	Math::Vector n = get_normal(Math::Vector());

	if ((std::abs(n.x) >= std::abs(n.y)) && (std::abs(n.x) >= std::abs(n.z))) {
    	alpha_11 = ((E_02.y * E_03.z) - (E_02.z * E_03.y)) / n.x;
    	beta_11  = ((E_01.y * E_02.z) - (E_01.z * E_02.y)) / n.x;
  	}
	else if ((std::abs(n.y) >= std::abs(n.x)) && (std::abs(n.y) >= std::abs(n.z))) {  
    	alpha_11 = ((E_02.z * E_03.x) - (E_02.x * E_03.z)) / n.y;
    	beta_11  = ((E_01.z * E_02.x) - (E_01.x * E_02.z)) / n.y;
  	}
	else {
    	alpha_11 = ((E_02.x * E_03.y) - (E_02.y * E_03.x)) / n.z;
    	beta_11  = ((E_01.x * E_02.y) - (E_01.y * E_02.x)) / n.z;
	}
	if (std::abs(alpha_11 - 1.0) < epsilon) {    

    	// Q is a trapezium.
		u = alpha;
		if (std::abs(beta_11 - 1.0) < epsilon) v = beta; // Q is a parallelogram.
		else v = beta / ((u * (beta_11 - 1.0)) + 1.0); // Q is a trapezium.
	}
	else if (std::abs(beta_11 - 1.0) < epsilon) {

    	// Q is a trapezium.
		v = beta;
		u = alpha / ((v * (alpha_11 - 1.0)) + 1.0);
	}
	else {
		double A = 1.0 - beta_11;
		double B = (alpha * (beta_11 - 1.0))
				- (beta * (alpha_11 - 1.0) - 1.0);
		double C = alpha;
		double D = (B * B) - (4.0 * A * C);
		double Q = -0.5 * (B + ((B < 0.0 ? -1.0 : 1.0)
				* std::sqrt(D)));
		u = Q / A;
		if ((u < 0.0) || (u > 1.0)) u = C / Q;
		v = beta / ((u * (beta_11 - 1.0)) + 1.0); 
	}	
	intersection_point = _ray.get_origin() + _ray.get_direction()*t;
	return make_tuple(infinity,intersection_point,this);

}

tuple<double,double> Quadrilateral::get_u_v(const Math::Vector& _point) const {
	double alfa = 0.0,beta = 0.0;
	
	double x = _point.x;
	double y = _point.y;
	double z = _point.z;
	
	double Xc = vertices[0].x;
	double Yc = vertices[0].y;
	double Zc = vertices[0].z;
	
	double Xa = (vertices[3] - vertices[0]).x;
	double Ya = (vertices[3] - vertices[0]).y;
	double Za = (vertices[3] - vertices[0]).z;
		
	double Xb = (vertices[1] - vertices[0]).x;
	double Yb = (vertices[1] - vertices[0]).y;
	double Zb = (vertices[1] - vertices[0]).z;
	
	alfa = (Xb*(y-Yc)+Xc*Yb-x*Yb)/(Xb*Ya-Xa*Yb);
	beta = -(Xa*(y-Yc)+Xc*Ya-x*Ya)/(Xb*Ya-Xa*Yb);
	
	return make_tuple(1-alfa,beta); /** FIX Me */
}

Quadrilateral& Quadrilateral::operator=(const Quadrilateral& _quadrilateral) {
	Object::operator=(_quadrilateral);
	vertices[0] = _quadrilateral.vertices[0];
	vertices[1] = _quadrilateral.vertices[1];
	vertices[2] = _quadrilateral.vertices[2];
	vertices[3] = _quadrilateral.vertices[3];
	return *this;
}

}

}
